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Recently we have introduced a simplified model of ecosystem assembly ( Capitan et al. 2009') for which we are able to map out all 
assembly pathways generated by external invasions in an exact manner In this paper we provide a deeper analysis of the model, 
obtaining analytical results and introducing some approximations which allow us to reconstruct the results of our previous work. 
In particular, we show that the population dynamics equations of a very general class of trophic -level structured food-web have 
an unique interior equilibrium point which is globally stable. We show analytically that communities found as end states of the 
assembly process are pyramidal and we find that the equilibrium abundance of any species at any trophic level is approximately 
inversely proportional to the number of species in that level. We also find that the per capita growth rate of a top predator invading a 
resident community is key to understand the appearance of complex end states reported in our previous work. The sign of these rates 
allows us to separate regions in the space of parameters where the end state is either a single community or a complex set containing 
more than one community. We have also built up analytical approximations to the time evolution of species abundances that allow 
us to determine, with high accuracy, the sequence of extinctions that an invasion may cause. Finally we apply this analysis to obtain 
the communities in the end states. To test the accuracy of the transition probability matrix generated by this analytical procedure 
for the end states, we have compared averages over those sets with those obtained from the graph derived by numerical integration 
of the Lotka-Volterra equations. The agreement is excellent. 
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1. Introduction 

A piece of common wisdom in ecology is that biodiver- 
sity enhances the stability of ecosystems. This has tradition- 
ally been a well established observational fact since the works 
of |Odum| ( fT953l ), |MacArthur| ( fT955| ) and [Eltonl ([1958} who 
showed that simple ecosystems (e.g. man-cultivated lands) un- 
dergo very large fluctuations in population and are vulnerable to 
invasion, an effect that gets reduced upon increasing the num- 
ber of predators and preys in the system. But early in the 70s 
May showed that randomly generated dynamical models for 
the populations of a community exhibit the opposite feature: 
the larger the species abundance the smaller its linear stability 
( [May 1972 1973 1. Thanks to this controversy we have gained 
very much insight into the nature of ecosystems ( [McCann] 
2000| l. Apart from the introduction of more refined concepts of 
ecosystem stability ('Pimm 1982 1, one of the main conclusions 



arising from the comparison of empirical data with May's pre- 
dictions on the bounds for community stability ( [Dunne 2006| l 
is that real ecosystem are within the tiny set of stable ones, no 
matter how large they are; in other words, ecosystems are far 
from being just random gatherings of species. 
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Natural communities carry out a selection mechanism that in- 
duces colonizers adaptation. There has been a lot of theoretical 
work in the past devoted to study the assembly of communities 



through successional invasions (Post and Pimm 
[19901 ICasel [T9901 |Law and Mo'rton[ |1993| 11996 



1983; Drake 



Morton and| 



Law 1997| l. Overall, these papers have provided a theoretical 



framework to understand how communities are built up ( iLaw[ 
1999| l. The basic process in which these models are based is 



the sequential arrival of rare species (invaders) that colonize 
the ecosystem and that may be established, possibly causing 
a global reconfiguration of the community in the long term by 
means of several species extinctions. Obviously, these mod- 
els are but idealizations of the complex processes taking place 
in real community assembly, but simple mechanisms acting in 
these models could be expected to be the ones responsible for 



the formation of real ecosystems ( |Law 1999 1. This approach of 
devising theoretical paradigms for real situations has been suc- 
cessfully applied over and over in the field of statistical mechan- 
ics — where, for instance, using such an idealization as the Ising 
model provides the clues to understanding ferromagnetism in 
real materials ( [Huang 1987 1. 



Previous assembly models tend all to rely on the Lotka- 



Volterra dynamics [but see the recent work of Lewis and Law 
( 2007| l], although differ in the criterion to accept an invasion. 



While Post and Pimm (1983i assumed that new species were 
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created ad hoc, according to certain stochastic rules, subsequent 



approaches (Drake 1990 Law and Morton 1996 1 introduced 



the concept of "species pool". A regional species pool is a set 
of possible invaders whose trophic interactions have been de- 
termined in advance ( [Law and Morton[ |1996| l. Despite these 
differences, all previous papers arrive at the conclusion that the 
species richness of each resident community increases along 
successional time, although the average resistance of a commu- 
nity to be colonized increases in time. Therefore community 
assembly increases biodiversity as well as stability, understood 
as resistance to invasions. 

Nevertheless, one must bear in mind that not all assembly 
pathways have been explored in these models. The conclusions 
reached so far rely on averages of quantities under study over 
a finite set of realizations of the underlying stochastic process, 
that is ultimately based on a finite pool of possible invaders. 
This has raised several question that remained without a defini- 
tive answer. For example, there was no clear-cut answer regard- 
ing the dependence of the results on the history of invasions. 
[Morton and Law| ( |1997[ ) found a final end state resistant to in- 
vasions by the remaining species in the pool at the end of the 
process, and this end state could be either a single ecosystem 
or a set involving more than one community connected by inva- 
sions with one another Despite this conclusion, the dependence 
of the end state on the assembly history is a matter of discussion 
( [Fukami and Morin] |2003| l. Moreover, we should not forget 
that the number of species in the pools employed is always rel- 
atively small, so the question remains as to whether larger pools 
lead to qualitatively different results. In this respect, it has been 



pointed out ( Case 1991 Levine and D' Antonio 1999 ) that the 
exhaustion of good invaders in the early assembly might be just 
an artifact of the finiteness of the pool. 

Trying to overcome the shortcomings of previous models, in 



our previous work Capitan et al. ( 2009^ ) we proposed a minimal- 
istic model of ecosystem assembly with which we were able 
to analyze all assembly pathways, thus characterizing the full 
assembly process. In spite of its simplicity, we recovered the 
same conclusions found previously. Our model is also based 
on a pool of species and a niche variable (the trophic level) that 
determines their interactions. In contrast, however, our pool is 
infinite. In spite of that, within the assumptions of the model, 
we found a finite number of (viable) communities linked by col- 
onization. This allowed us to define an assembly graph for our 
model — similar to that of [Warren et al.| ( f2003| , who studied 
the assembly process experimentally for a small pool of 6 pro- 
tist species. By assigning transition probabilities to the links 
of this graph the assembly process was mapped to a Markov 
chain ( [Karlin and Taylor| [1975| l, which is tantamount to say- 
ing that we defined a statistical mechanics on the set of viable 
communities (microstates). In other words, our model gives 
the probability distribution of all these microstates at any time. 
This allowed us to characterize both transient and equilibrium 
states, as well as to compute the time evolution of any observ- 
able along the assembly in an exact manner. But more impor- 
tantly, as our model provides a complete and exact (albeit nu- 
meric) description of the assembly process, we can positively 
state that, under the assumptions of our model, in the long-term 



assembly dynamics a unique enstate is reached, and this state 
is formed by just one uninvadable community or a closed set of 
communities connected between them. These sets contain the 
communities that survive in the long term, and the ecosystem 
can be regarded as a fluctuating community that can vary each 
level occupancy trough successional invasions. 

In this paper we will give some analytical results for the un- 
derlying population dynamics of our assembly model, and we 
will see how these results can be combined together to arrive at 
the same conclusions we obtained numerically in our previous 
work. Relying on these analytic results, we will be able to de- 
scribe the observables that characterize the end states with high 
accuracy. In particular, we will reproduce the variation of the 
number of communities in each end state with the abundance 
of abiotic resources, as well as the average values of quanti- 
ties like the species richness. We will leave the computational 
and numerical results that can be obtained with this model for 
the second paper of this suite (Capitan et al. 2010), which will 
be focused in the successional variation of biologically rele- 
vant quantities along the assembly, and the analysis of the main 
properties of transient states. 

This paper is organized as follows. Section [2] is devoted to 
the analysis of a rather general model of trophic-level struc- 
tured food-webs, and the discussion of its dynamic stability. In 
Section[3]we will restrict ourselves to a particular case of com- 
munity by making a species symmetry assumption, that renders 
our model closer to neutral models and allows a more detailed 
analytical study. In Section [4] we will deduce some analyti- 
cal properties of the equilibrium point, such as estimations of 
the maximum number of species allowed in a community for a 
given set of parameters, or the maximum number of trophic lev- 
els that the amount of resource allows. Section|5]is dedicated to 
discussing some criteria for an invader to establish in a commu- 
nity, and to give some global analytical approximations to the 
time evolution of a system invaded by a top predator. Finally, 
in Section [6] we will apply our analysis to recover the results 
obtained in Capitan et al. (2009) by means of a numerical inte- 
gration of the population dynamics equations. 

The two papers of this suite are self-contained and can be 
read separately, although they are cross-referenced. Readers 
interested in the underlying population dynamics of our model 
will find a detailed discussion in this paper Those readers more 
interested in the ecological consequences and results that the 
model provides can skip the technical Sections [4] and [5] For a 
full account of the results that we have obtained, we refer them 
to the companion paper. 

2. Trophic-level structured food-webs 

How species are arranged in a network to conform a food- 
web is a question difficult to answer. The specific topology 
of the network where feeding interactions take place is very 
complex and several complicated models have been proposed 



for both the structure and the dynamics of food- webs ( Dunne 
2006). In contrast, our aim in Capitan et al. (20091 was to 



construct a minimalistic model, so we considered the tradi- 
tional picture of trophic pyramids of interacting species in dif- 
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ferent, well defined trophic levels. Although trophic levels can 
be roughly described in real webs ( Martinez et al. 2006| l, we 
will assume that feeding interactions take place strictly between 
species belonging to contiguous, well defined trophic levels. 
This is a standard (and accurate) assumption, as the models of 



tri-trophic food chains show (Bascompte and MeUan 2005 i. 
This notwithstanding, it is acknowledged that omnivory, i.e. 
predation from several levels, exists although is still an open 
question how common it is. For example, work on food-web 
motifs has found that omnivory is sometimes under-represented 
and sometimes over-represented in real networks ( |Bascompte| 
and Melian 2005 | l. However, the impact of including omnivory 
in the model could lead to non trivial results. Since the trophic 
level is normally related to species size, feeding from lower lev- 
els will provide less energy to predators, so proper allometric 
relations should be included in the model to fix the interaction 
strengths. For the sake of simplicity, we will not divert our- 
selves from the standard assumption of disregarding omnivory. 

Therefore, any species at level £ will feed only on species at 
level i - 1 and will be predated only by species at level £ + 1. 
Let sc be the number of species in the ^-th level. Thus for an 
ecological community with L trophic levels the total number of 
species is 5 = Yjt=i-^{- I" order to determine which species are 
predated at each level, we define the set of interaction matrices 



the species at trophic level £. Following BastoUa et al. ( 2005a i 
we propose the mean-field dynamics 



natural way to represent this matrix is 



(2) 



where < 1 measures the relative magnitude between intra- 
and interspecific competition, and I is the identity matrix. Di- 
agonal elements of are equal to 1 due to the normalization of 
the intraspecific competition. We will assume (the reasons will 
become clearer later) that the competition matrix is symmetric 
and positive definite. 

Interspecific competition due to sharing common preys is im- 
plicitly represented in the predation terms. There is however a 
direct competition due to other effects, such as territorial com- 
petition, mutual aggressions, etc. We will assume [as in Bas- 
itolla et al., (2005b) ] that species sharing more preys are closely 
related ecologically [this fact might have support from a evolu- 
tionary viewpoint as shown in Rezende et al. ( 2007| l], so their 
requirements are similar and we can assume that elements of 
are proportional to the ecological overlappi ng between species 
(iLassig et al.| |200ll iBastoUa et al.| |2005b') 



F , with dimensions S( x Sf_i, such that the element F^ - 1 so that 
when species j in level - 1 is a prey of species / in level £, and 
is zero otherwise. Any particular choice of this set of matrices 
determines the food-web in our model. 

According to our aim of developing a simplified model, we 
propose a simple population dynamics with the purpose of cap- 
turing on average the main behavior of species abundances. It is 
inspired in a model used before to study coexistence in compet- 
ing communities ( [Lassig et al.||200l]|Bastolla et al.||2005a|b| l. 
Population dynamics is modelled by Lotka-Volterra equations, 
including both predator-prey interactions as well as intra- and 
interspecific competition. Thus, in order to keep the model min- 
imalistic we have chosen not to include other interaction types 
such as mutualism. 

Let be a column vector with the population densities of all 



Let jT^.j represent 
the number of corntnon preys tor species / and j belonging 
to level £. The species overlapping due to common preys is 

K^j - n^.jl ^7r[;7rj, with 7r[ the total number of preys of species 
/. Under our matrix notation, n^.. - (F^F'^^),-,- and n^. - (F^F^^),/, 



= (1 - pf)l + p^dV(dV)'^ 



(3) 



being a diagonal matrix with elements (F'^F^ . Ex- 
pressed as ([3|, it is evident that such a competition matrix is 
symmetric and positive definite. It is worth mentioning that 
this system does not fulfil the hypotheses leading to Cause's 
competitive exclusion principle ( [Hofbauer and Sigmundj 11998} 
IBastoUa et al.||2003al l, even when there is a single level. Among 
other things, this is due to the fact that competition coefficients 
between different species are not all the same. This point will 
be discussed in more detail in the second paper of this suite 
( |Capit anet"aLl [20T0l l. 

We regard all the species as consumers, and so they have a 
death rate, a^., which is the i-th component of vector a^. Note 
that in a real food-web the interaction coefficients will not be 
uniform within a trophic level. In this sense, we represent inter- 
actions averaged (mean-field) in each level but we allow varia- 
tion in the strength of the interactions among different trophic 
levels. Finally, all species at the first level predate on a single 
resource, whose time evolution is given by 



We assume that the strength of the feeding interactions between 
contiguous levels is fixed and determined by the constants y'^, 
which control the amount of energy available to reproduction 
for each predation event for species at level {, and y'^ (> y+), 
which take into account the mean damage caused by predation 
over level {. The ratio y+ly''_ measures the efficiency of conver- 
sion of prey biomass into predator biomass. 

Interspecific competition in a trophic level is measured by the 
off-diagonal elements of the si x S( matrix , while intraspe- 
cific competition (diagonal elements) is normalized to unity 
(this just amounts to fixing a time scale for the dynamics). A 



(4) 



The constant R is the maximum amount of resource in the ab- 
sence of its consumers. The abundance rP has to be understood 
as the amount of a primary abiotic resource, like sunlight, wa- 
ter, nitrogen, etc. It has to be considered as an energetic input 
for the maintenance of the remaining species in the community. 
The amount of such resource is limited, hence the saturation of 
rP at a value R. 

The model is supplemented by an extinction threshold, «e > 
0, independent of the species. If a population falls below this 
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value it is considered extinct (real populations can not be arbi- 
trarily small). This viability condition has been previously used 
in similar models ([Kokkoris et~aLj [T9991 |Borrvall et aT] [20001 
[Eklof and Ebenman, ,2006^ , and accounts for the vulnerabil- 
ity of low density communities against external environmental 
variations or adverse mutations (Pimrnj |1991| l. The technical 
need for this extinction threshold in our model will become 
clearer when we describe the variation of the densities in terms 
of the occupancy of each level. 

2.1. Dynamic stability of the interior equilibrium point 

Equations ([TJ, Q have several equilibria. Among them, the 
main one is obtained by equating the right-hand side of these 
equations to zero. If all the equilibrium densities are positive, 
this fixed point is called the interior equilibrium. The popula- 
tion at equilibrium are obtained as the solution of the linear 
system of 5 -H 1 equations 



yir'p'-' - B'p' - y'_-\T'^'fp'^' = a\ 
p'^ + yKrVp'^R. 



(5) 



for £ - 1, . . . , L. The remaining equilibria are obtained by set- 
ting to zero any subset of the populations and solving the Unear 
system resulting from eliminating those variables. The result- 
ing system is the same as (|5]l but if species / at level { has zero 
equilibrium abundance, the i-th column in the corresponding 
matrix has to be eliminated. Therefore one only needs the 
solutions of the linear systems Q for a given choice of the set 
of matrices {T'^j^^j in order to fully determine all the equilib- 
rium densities. 

Since feeding relations are established among contiguous 
levels, (|5| acquires a block-tridiagonal structure. Due to this 
form, the interior equilibrium can be formally obtained by ap- 
plying Gaussian elimination. We put the equilibrium abun- 
dances in the form 



(6) 



for certain x sc matrices and Sf^i x 1 vectors to be 
determined (£ = 1,. . . ,L + 1). Substitution into Q gives the 
following recursive relations for and c^. 



(7) 



Since the resource can only be predated and there is no compe- 
tition, we set r'^ = and p" - 0. This leads to the initial condi- 
tions M' - -7'i(r')^ and c' = -7? according to Q. Thus, given 
a particular set of matrices {T'^l^^j, (|7]i fully determines and 
c^. After that, starting from the boundary condition p^^' = 
(the community has exactly L trophic levels), we backsubstitute 
in (|6| to get the equilibrium densities. 

We can push further the property that our dynamical sys- 
tem ([T| is block-tridiagonal to study its dynamic stability. Let 
us show that interior equilibria p^., for all / = 1 , . . . , if and 
£ = 0, . . . , L, are globally stable. This result is based in the 
existence a Lyapunov function ( Hofbauer and Sigmund 1998[ l, 



which guarantees that any positive initial condition evolves to- 
wards the interior equilibrium. The Lyapunov function for this 
system is 



(8) 



where Ai 



=1 



^ for ;t 

n 



1, 



, L and A 







1. 



For (|8]l to be a Lyapunov function, we just need to check that 
•y < along any orbit {«^(0}^=o starting with positive initial 
abundances ("Hofbauer and Sigmund 1998 1. Let us compute its 



time derivative. If we consider the displaced variables 



we can write ([T]l as = "f^f' where 

^yir^-' -By-yi{r'^'f 



(9) 



(10) 



hence the time derivative is simply 'V({«'^)) 
Yj{=o^e y] After substituting ([T0|, we arrive at 



^({n'))^-Y^At{y'fE 



{=0 
L-1 



(11) 



Thus our previous choice of A^ cancels the second sum. Since 
is positive definite, we deduce that the time derivative of the 
Lyapunov function is negative along any orbit, and therefore 



Lyapunov's theorem (Hofbauer and Sigmund 1998) ensures 
the global stability of the non-trivial rest point p*-. Note that 
the existence of this Lyapunov function is a direct consequence 
of the block-tridiagonal structure of the dynamical system ([T])- 
(|4]|, hence the assumption of predation only between contiguous 
levels ensures this global stability. 

3. Species symmetry assumption 

In what follows, we will restrict ourselves to the dynami- 
cal system ([TJ with the particular choice of interaction matrices 
F^. - 1 for any /, j, (. This was the system studied in Capitan 



et al. (2009). This assumption implies that all species are gen- 



eralist, and the model can now be regarded as a mean-field-like 
picture of real communities, since all species in contiguous lev- 
els interact with each other We will assume as well that in- 
teraction coefficients are independent of the trophic level, and 
we will simply denote them as y+, y^, p and a. These pa- 
rameters should now be understood as an average strength of 
the processes involved in the population dynamics. These kind 
of models, which do not make any explicit difference among 



species, are referred to as neutral (Hubbell 2001 Etienne and 



Alonso 2007 1. From the point of view of the trophic interac- 



tions there is no difference between species (neither the rates 
nor the set of preys they feed on make any distinction among 



species). We introduce this symmetric scenario because it will 
allow a simpler, analytical description of the community. 

Pure neutral models do not make any distinction whatso- 
ever between species. This is not our case, because species 
can be distinguished by their different balance between intra- 
and interspecific competition. Neutrality in our model has to be 



understood as a species symmetry assumption (Alonso et al. 



2008 1 for the strength of the interactions. We will discuss the 
case p - I, when the model turns to be fuUy symmetric (strictly 
neutral), in the second paper of this suite ( [Capitan et al. 2010 1. 

Under this symmetry assumption, the population dynamics 
(|T|l with the competition matrix (|3]l transforms into n[ = ^fnf. 
where 



-a + y+N^-^ - (1 - p)nf - pN'' - J-N' 



1i 



(12) 



being A^'^ = YfiL\ ■ The set of equations Q for the interior rest 
point imply that the equilibrium abundances are equal for any 
two species / and j of the same level. Hence the equilibrium 



abundances {p^]^f, are the solution to the linear system 



a = y+se-ip - [1 +pisc-l)]p -y-SMp'^\ 
R^p"+y-sip\ 



(13) 



for ( - I, . . . ,L. Note that the global stability result holds only 
for this equilibrium point. 

3.1. Reduced dynamical system 

As in our previous work (Capitan et al. 2009'), equilibrium 
communities will undergo invasions. Thus we are interested in 
the time dynamics of an invaded community initially at equi- 
librium. Notice that the per capita growth rates ( [T2] i satisfy the 
equality 



.,n,.. 



under the interchange of the abundance of two species at the 
same level. This symmetry, together with an initial condition 
where nf(0) - n^(0), is enough to show that the time evolu- 
tion of both species is identical (see [Appendix A| i. Thus we 
can reduce our dynamical system to a set of L H- 1 differential 
equations. 



= -a + y+se^in'^ ' - [1 + p(se - l)]n' - y^se+in'^\ 



^ R-n^ - y^sin. 



(15) 



There is another crucial difference between our model and 
usual neutral models in the literature. Although neutral models 
ignore species identity, they are stochastic. It is the ecologi- 
cal drift what makes species abundances to stochastically vary. 
This stochasticity is the ultimate reason for extinction in neutral 
models. On the contrary, our dynamical system is determinis- 
tic. The reason to include the (somehow arbitrary) extinction 
threshold is to "mimic" this fluctuation-driven extinction of 
species with low abundance. 



Thus extinctions must be understood stochastically in our 
model. As it was pointed out in Capitan et al. ( 2009[ l, the 
stochastic effect of adverse mutations or external variations of 
the environment that make species to go extinct is taken into 
account in our deterministic dynamics with the viability con- 
dition > ric. Notice however that, strictly speaking, when 
a species of one level falls below the whole level does too. 
Extinguishing the whole level as the strict dynamics would re- 
quire would be unrealistic. Instead we eliminate species one 
by one until viability is recovered (Capitan et al. 2009]l. This 
latter dynamics would approximate better what one would find 
in a truly stochastic neutral model, in which the simultaneous 
extinction of general species is very unlikely to happen. 

3.2. Structural stability 

We have chosen the constants to be uniform in our model, 
this making all species on each trophic level at equilibrium have 
equal abundance. However, according to competitive exclusion 
(MacArthur and Levins 1964| l, a tiny variation in the param- 
eters that makes any difference among species will make the 
system unstable. Fortunately, for this class of models the com- 
petitive exclusion principle does not hold as such. This has been 
discussed at length in |Bastolla et al.| ( [2005a[ ). In this paper the 
authors derive some bounds to the variation allowed for the con- 
stants that the system can tolerate without leading any species 
to extinction. In fact, the dynamical system they discuss is the 
same as we have described in Section |2] with different con- 
stants for different species. The more diverse the ecosystem is 
the stricter are these bounds, but in any case, no matter how 
diverse the ecosystem is, some variation of the constants is al- 
ways tolerated without this leading any species to extinction. 
This proves the structural stability of our system, even under 
the assumption of species symmetry. 



) (24) 4. Analytical properties of the interior rest point 



4.1. Maximum number of species and maximum number of lev- 
els 

In this subsection we will obtain an analytical estimation of 
the maximum number of species that a trophic level can host 
among all the possible viable equilibria. We simply set all the 
abundances in each level to be equal to «c and solve the result- 
ing linear system ( [T3] l for {sc]^_^ and S(, = p^ jn^. 



R 



So + y-s\ = 



y+if_i - psc - y-se+\ = 1 - p -t- 



(16) 



for { > 1. We introduce the generating function G{z) = 
Yi7=()^(^^ for the sequence {si]^^y The explicit solution will 
depend on two initial conditions so and si, since we have a 
two-term recursion. We will leave them undetermined for the 
moment. The second equation of ([T6| allows us to calculate 
explicitly G(z), 



G(z) = 



(1 -p + a/nc)z^ 



y-So+z(psQ +y-si) 



(1 - zXy+z^ - pz-y-) y+z^ - pz-y- 



■ (17) 
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We recover the general term of S[ by a series expansion of the 
generating function. Let us first define the constants fi = (1 -p+ 
o^/nc)/(7--y++p) andzi = (p+ + 4y+7_)/(2y+). In order 
to get compact expressions, we define the auxiliary sequence 



a{ = 



z+ - Z- 



(18) 



which satisfies the two-term recursion j-Uc - pai^x + 7+af-2 
with initial conditions fl_i = 0, oq = 1- This recurrence can be 
fully expressed as a linear combination of powers of p/y_ and 



at 



Lf/2J 

z 

/t=0 



(19) 



for all ^ > 0, [jcj denoting the integer part of x. 
Expanding G(z) we obtain st in terms of at. 



(-1)' 



r- 



(io + yu)af-2 - (*1 + /^)flf- 



(20) 



for £ > 2, where a( can be evaluated either using ([TSj or ([T9|. 
In order to solve the system ([T6]l, we have to impose ,s-£+i = 
for an ecosystem to have L trophic levels. This provides a 
linear relation between and s\ which, together with the first 
equation of ( [T6] l, forms a linear system that determines both sq 
and s\. The result is 



{Rjuc + fxj- + n)aL - (- 1) Vr- 
*o = 

ai + r+fli-i 

y+(Rln,+ ny_+ ^j)aL-i + (-1)^7- 
s\ = z — ^ Z^- 

Substituting ( |2T| ) into (|20| and taking into account that 



(21) 



flZ,flf_2 - flL-lflf-I = (-1) 



r- 



t-i 



ai-t 



is a direct consequence of the recurrence satisfied by 
finally get 



ai-t 



ai + y+ai-i 



(22) 



we 



(23) 



for all ^ > 1. This is the analytic solution of the system ([T6| 
and gives an estimate of the maximum occupancy per level as a 
function of the parameters of the model. Note that, despite what 
( fTSj ) might suggest, no additional factors of the form y+ly- can 
be extracted from ac according to ([T9]l, so the lowest power of 
the ratio y^ly- in the expression for S[ is (y+Zy-Y- 

This dependence of se on (y+Zy^Y is remarkable. In fact, 
in our previous work (Capitan et al. 20091 we observed that 
the communities in the end states of the assembly process were 
pyramidal. This is, in turn, a consequence of the exhaustion of 
the species occupancy in each trophic level. Notice also that the 
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Figure 1: Approximate equilibrium densities. Starting from a community with 
4 levels and occupancies .vi = 127, S2 = 58, S3 = 7 and S4 = 1, we plot the 
variation of as a function of if, which exhibits a dependence C/.Vf . Full lines 
with circles show the exact solution of l |13) , and dotted lines with crosses show 
our approximation \2(s\ . Insets contain the relative error of each approximation. 
Remaining parameters are R = 1505, y+ = 0.5, 7_ = 5, p = 0.3 and a = \. 



estimation of the maximum number of species that a commu- 
nity can host depends linearly on the resource saturation. This 
linear dependence on R was also observed in our previous work. 

Our estimation of the maximum occupancy of each trophic 
level also provides a condition for the maximum number of 
trophic levels that a set of parameters allows. Imposing si> \ 
yields a condition for the allowance of L trophic levels. 



R jy 

-+Mr- + i)> - 
"c \y 



[(1 + ^i)iaL + y+ai-i) 
+ IJ{aL-\ + y+ai-i)] ■ 



(24) 



Therefore we have a minimum value of the resource saturation 
for L trophic levels to be viable in a community. 

4.2. Approximation of the equilibrium abundances 

In our model, each set {sc]^^q of species occupancies in each 
level determines a set of equilibrium densities according to ( [T3| l. 
Finding p^{{sk]) is difficult, but in this section we will give a 
rather good approximation for large enough S(. First we write 
the system in terms of the total population at each level, - 
Sep' (^= 1,...,L), 



(25) 



+ y_P^ = R. 



Written in this way, it seems natural to expand the solution in 



powers of 1 1 s. In Appendix B we show that we can approxi 
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mate 



As we can see in Figure [T] this first order approximation cap- 
tures accurately the variation of the equiUbrium densities p'^ 
with S(. Besides, we also obtain a very accurate approxima- 
tion when we vary the number of species sj in levels other than 
£. Note that, even when the occupancy of a level is small (lower 
panels of Figure [TJ, the approximation remains good. 

In the limit S( » 1 we obtain the dependence x C/sg, 
which provides the general tendency observed in Figure [T] 
Moreover, in the biologically relevant limit R » a, and tak- 
ing into account the explicit expressions for Tij and given 
inT 



(26) 



Appendix B populations behave like 



R Ir 



ai-t 



s{ \y-i flL + r+flL-i 



(27) 



for ( > 0. Several conclusions can be extracted from this de- 
pendence. First, when the number of species in the ^'-th level 
is exhausted, according to ( |23] l, we obtain a population density 
p^ ^ He, as expected. But more importantly, it represents an- 
other reason for the extinction threshold to be included in our 
model. If there were no threshold, equilibrium densities would 
monotonically decrease with without ever becoming zero. 
The assembly graph would then contain infinitely many com- 
munities thus becoming intractable. 

5. Invaded dynamics 



In Capitan et al. ( 2009 1 it was assumed that, during the as- 



sembly process, successional invasions occur and modify resi- 
dent communities at equilibrium. There we made the hypoth- 
esis of the average time between consecutive invasions being 
much longer than the typical dynamic time scale for the com- 
munity to reach the equilibrium state. This is actually what is 
observed. In relation to the different time scales between in- 
vasion and competition, invasion events may take place at the 
scale of decades, long enough time for invaded communities 
to stabilize [for example, the rate of new invasions in islands 
may be one every few year (Sax et al. 2005| l]. This assump- 
tion has also been made in previous papers Uke Kokkoris et al. 
( |1999| l, where authors assume that after each invasion there is 
a re-organization of the community prior to a new invasion. 
Specifically, they solve the dynamical system describing the 
new community with the invader until reaching the carrying ca- 
pacity. These new densities are then used as initial values for 
the new systems resulting from the next invasion [see details 
in [Kokkoris et al.| ( |1999| )]. The same idea was applied in the 
construction of our assembly model (Capitan et al. 2009[ l. 

We used a second hypothesis as well, namely that the popu- 
lation of the invader is small (equal to the extinction threshold 
ric). This is what is actually found in real situations. It is a 
well established fact that colonizers rarely reach a new habitat 
in high numbers ( Roughgarden 1974 Turelli 1981[ l. In theory, 
the probability of a small propagule to extend is used as the 



invasibility criterion. In biological control, management of in- 
vasions is based on looking for a small density of species in new 



areas (Liebhold and Bascompte 2003 i. In this case, theoretical 



and empirical work has taken advantage to predict conditions 
of eradication based on density thresholds (Allee effects) and 
demographic stochasticity. 

Therefore we can assume invaders arriving at some level of 
a community in equilibrium with a small abundance set equal 
to the extinction threshold. Under the species symmetry as- 
sumption, the dynamic system h^. - n^.q^j given by the re- 
sponse function ( [T2| ) applies as well for the invaded system, 
with A^'^ = YjILi nf+n and n being the population density of the 
invader. Therefore, once the equilibrium is reached after the in- 
vasion, the density of the invader will equal p'- (the density of 
the remaining species in that level), which can be obtained by 
solving ^T5\ with an occupancy -H 1 in the ^-th level. More- 
over, the global stability condition applies as well to the invaded 
dynamics. So we just need to check the viability of the resulting 
equilibria in order to determine whether the invader is accepted. 

If the invasion takes place at level L + 1, the equation for the 
invader is simply 



- — —a + j+siti - n, 
n 



(28) 



which in fact is the last equation of the system ( [T5| ) for a com- 
munity of L H- 1 levels with occupancies {si, . . . ,sl, 1). Hence 
the global stability condition still remains applicable and the 
invader will be accepted if the resulting equilibrium is viable. 

The complexity of the assembly dynamics comes from the 
cases where some level in the invaded community falls below 



the extinction threshold. The approach we used in Capitan et al. 



( 2009 1 to determine the sequence in which species go extinct 
until leading to a final viable ecosystem was the following: for 
the levels that fell below the extinction threshold once the equi- 
librium had been reached, we went back in their trajectory to the 
point where the population of some species crossed the extinc- 
tion level He for the first time, we eliminated one species from 
that level and restarted the dynamics from that point. In this 
paper we will propose an alternative way to determine that se- 
quence based on several criteria and analytical approximations 
that we will discuss below. 

5.7. Invasion criteria 

Consider the general dynamical system ij/jc,- - qi(x,xj), 
i//x/ - qj(x,xi) for an arbitrary community with S species, 
where x are the densities of the species in the resident commu- 
nity and xj is the density of the invader. The establishment of a 
colonizer in systems of this kind depends crucially on the initial 
per-capita growth rate of the invader (Law and Morton 1996| l. 
In fact, the condition that must be satisfied for a new species to 
increase when rare is 



lim 



r Jo 



qi(x(t),Xi - 0)dt > 0, 



(29) 



i.e., the time average of the per-capita rate of increase of the 
invader is positive when the species of the resident commu- 
nity remain under certain attractor x(t) of the dynamics. In our 
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model, the only attractor is the interior rest point, so the con- 
dition reduces to qi{p, 0) > 0, where p is the rest point of the 
resident community. Strictly speaking, our model has a non- 
zero extinction threshold, so this condition has to be replaced 
by qiip, ric) > 0. Since we start from a resident community 
initially at equilibrium and the invader initial density is n^, this 
condition reduces to the initial per-capita growth rate of the in- 
vader 

The condition qi(p,nc) > can be used to obtain criteria for 
the invasibility at each level. For example, consider the initial 
growth rate of the invader when the invasion takes place at the 
level L + I [Eq. (|28|l]. The condition for this rate to be positive 
is 

ry 4- n„ 

(30) 



If this condition does not hold, the invader is the first species 
to go extinct because it starts at the extinction level and with 
a negative initial rate. In the end states, the populations of the 



resident community are close to (but above) «c (Capitan et al. 



2009 1, so the former condition provides the approximate bound 



a + Hr 



(31) 



Even if the initial growth rate of the invader is positive, asymp- 
totically the level L + \ may not be viable. If this happens, 
during the time when the population of the invader is above 
«f, extinctions may occur at lower levels. This situation ex- 
plains the accumulation of recurrent states that we observed in 



Capitan et al.| ( 2009 1 when we varied the resource saturation 
(see Section [6]l. 

Invasions at levels ^ < L are subject to similar conditions. 
For the initial growth rate of the invader to be positive 



1 



(32) 



must hold. In general, an initially positive growth rate could 
lead to potential extinctions in the remaining levels while the 
equilibrium density of the invader is above the threshold. But 
it could happen as well that the invader extinguishes at equilib- 
rium with some initial transient time above the extinction. To 
estimate a condition for this to happen, let us assume that den- 
sities and occupancies are inversely proportional (see (|26| and 
Figure [T]i. Then the equilibrium abundance of the invader is 
S[p^ l{_S( + 1), therefore if 



p <nAl + 



(33) 



the invader goes extinct. This condition, together with ( |32| l, 
leads to 

se<--\ (34) 
P 

so below this bound, the invader initially grows but becomes 
extinct at equilibrium. We wiU use this condition to explain the 
appearance of some recurrent subsets for certain values of R 
(see Section[6]l. 

It would be nice, however, to have a systematic way to pre- 
dict the sequence of extinctions after an invasion has occurred. 



Based on our approximations for the equilibrium densities, we 
can propose a way to sequentially remove species for invasions 
at lower levels. Within the end states of our model, abundances 
are close to the extinction threshold. Then ( |23| ) implies that 
communities are pyramidal, so lower levels are highly occu- 
pied but higher levels contain a small number of species. Ac- 
cordingly, the increase of one species in a lower level has no 
significant effect in the equilibrium abundances of the commu- 
nity. Therefore if a species goes extinct after an invasion in a 
low level, it has to be the invader itself. 

The extinction sequence for invasions in higher levels is not 
so easy to predict. Nevertheless, changes in abundances upon 
increasing st are larger the higher the level (Figure[T]) so, in case 
that several levels fall below the threshold, we can make the 
assumption that it is always the "highest" species the one that 
goes extinct first. This procedure provides a certain sequence 
of extinctions whose accuracy will be checked in Section|6] 

The prediction of the sequence of extinctions can be com- 
plicated when a top predator invades if the resource saturation 
values do not allow for L -H 1 levels. We have devised global 
approximations to the dynamics in this case to predict the order 
of extinctions without having to resort to the numerical integra- 



tion of the system of differential equations, as we did in Capitan 
|etaL] ( |2009] l. 

5.2. Global approximations to the dynamics invaded by a top 
predator 

Our heuristic approximations to the time dynamics of the 
system ([TSj when an invader arrives at level L -i- 1 are some- 
how inspired in the matching technique used to obtain analytic 
approximations to perturbed differential equations [see, for ex- 
ample Bender and Orszag ( 1984| l]. First we calculate the equi- 
librium point {p'^lf^g by either solving ( [T3| ) or using the approx- 
imations (|26|. Then we approximate n^+'(f) by the sum of its 
long-term dependence n^^'(0 (near equilibrium) plus a short- 
term behavior n^(^'(f). For the long term, a linear stability anal- 
ysis shows that the solution exponentially decays towards the 
equilibrium point, so we will set 

„^+i(f) = + e-^'[di) cos(wO + di sin(wO] (35) 

where the eigenvalue of the linear stability matrix which is clos- 
est to zero is -A + ico (cj may be zero). The constants t/o and di 
remain undetermined for the moment. 
For the short-term behavior we propose 



'(f) = C(t)e-^', 



(36) 



where C(f) - YjJ Cjt^ is a polynomial whose coefficients and 
the exponent ^ need to be determined to capture the transient 
time evolution. This way to express the short-term behavior is 
inspired in the initial transient decay that can be observed in the 
initial dynamics prior to getting close to the equilibrium point 
(see Figures [2] and |3]l. The polynomial has been included so 
as to properly capture the initial condition and the initial devia- 
tions to the exponential decay. The technical details to calculate 
the undetermined coefficients in ( (35| ) and ([36| are deferred to 
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Figure 2: Dashed lines show our approximation for the dynamics of a four- 
level community determined by the occupancies S[ = 1 10, S2 = 50, = 6 and 
i4 = 5 when invaded by a top predator at level 5. For this case the eigenvalue 
closest to zero of the linear stability matrix is complex. Full lines represent the 
numerical integration of \12) . Remaining parameters are the same as in Figure 
[T] The whole time evolution is accurately predicted. The extinction level = I 
is showed as a dotted line. We can see how the first extinction in the community 
takes place at level 4. 



Appendix C Figures |2] and [3] illustrate the validity of this ap- 
proximation in capturing the global trend of the time evolution. 

To reproduce the ordering of the extinctions we need the ex- 
tinction times for each level, and these times are approximated 
with a higher accuracy than the dynamic trajectories themselves 
(see Figure [3]). In Figure [4] we illustrate, for a particular com- 
munity, the extinction procedure compared to our analytical ap- 
proximations. In this case, the first level falling below rtc is the 
fourth one (upper panel). Then we remove one species from 
that level and restart the dynamics from the point of extinction, 
and the fourth level falls again below rtc (second panel). After 
the removal of a new species, the fourth level ends up above ric 
at equilibrium. Now the next level ending below «c is the sec- 
ond one. We move to the point of extinction of this second level 
and restart the dynamics after removing one species from £ = 2. 
After that it is just the invader ({ - 5) the only one that falls be- 
low the threshold, so we remove it and the resulting community 
becomes viable. Were it not, we would apply the same extinc- 
tion procedure again and again until the final community is vi- 
able. The sequence of extinctions is well reproduced with our 
approximate solution, although slight differences that alter the 
order of extinctions may occur when different levels fall below 
tic roughly at the same time. 



6. Application to community assembly 

Our goal in this paper was to provide analytical support. 



albeit approximate, to the results obtained in Capitan et al. 
( 2009[ l. We want to check now whether our approximations 



Figure 3: Same as Figure |2] but with R = 1200 and occupancies S[ = 106, 
S2 = 49, S3 = 6 and S4 = 4. For this case the eigenvalue closest to zero of 
the linear stability matrix is real. Although there is some discrepancy in our 
approximations, the global trend is captured and the extinction times after the 
invasion are accurately predicted. 
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Table 1: Estimation of the value of Rjric for the appearance of a recun'ent 
set with more than one community (left). Minimum values of that allow 
a community with L levels, according to j24j (right). The interval of values 
of R that correspond to the recurrent sets is approximately [Sicc,Smi„]. R\^^ 
and R' are the corresponding values found using numeiical analysis |Capitan| 
[erari(2'009/ mapping the whole range of R with a resolution AS = 5. 



assembly process. With this aim, we have varied the parameter 
R within the range from 10 to 1700 in steps A/? = 5. The re- 
maining parameters of the model will be set as in our previous 
work: y+ = 0.5, j- - 5, p - 0.3, a - \ and ric -\. 

Let us first fix the number of levels L. We can determine with 
(|24li the minimum value R.^,„ that allows L+ \ levels. The results 



are summarized in Table [T] Moreover, we can combine ( |23] l 
and ( |3T| ) to give an estimation of the initial value of R^^^ for the 
appearance of a recurrent set with more than one community. 



R /y_ 

-+Mr- + i)> - 



a + llr 



y+ric 



(37) 



correctly predict the recurrent sets which are end states of the 



The resulting values show a good agreement with those ob- 
tained numerically in Capitan et al. ( 2009) (see Table[T]l. 

Then, for a given R, we can read off from Table [T]the number 
of levels for the communities within the recurrent set. Once we 
know it, we determine with (|23ll an estimation for the maximum 



9 



1.5 \ 



1.2 
1 

0.8 
1.2 
1 

0.8 
1.2 



0.8 



level 2 
level 4 
level 5 



1000 




1 - 



0.5 



1.5 



Figure 4: Extinctions sequence for the community with si = 110, S2 = 51, 
S3 = 6 and i4 = 5 invaded at level 5 (parameter values are the same as in Figure 
[T] and «(. = 1 is showed with a horizontal dotted line). We just show the time 
evolution of the levels that go extinct or are close to extinction in equilibrium. 
Dotted curves con'espond to our analytical approximations. We show, with 
vertical lines, the time of the first level that go extinct. The sequence of extinct 
levels is 4, 4, 2, 5 until viability is recovered. 



occupancies allowed. We round off the estimates to get an inte- 
ger set of values {s(} and calculate the associated interior equi- 
librium. It can happen that some of the fall below tic, so we 
decrease the corresponding occupancies sc eliminating species 
one by one until the equilibrium turns out to be viable. This 
way we obtain a community very close to those of the recurrent 
set (communities within this set are close to extinction), so we 
can use it as the initial community to start the assembly process. 
We then compute the set of viable communities connected to it, 
which defines an assembly graph much smaller than those ob- 
tained in Capitan et al. ( 2009[ l starting from the empty commu- 
nity 0. We analyze the graph to obtain its recurrent sets using 
the algorithm of'Xie an d Beerel| ( |1998[ ) and we get one single 
set. In Figures |5] and [6] we plot the number of communities in 
each end state, showing a good agreement between the results 
obtained with the analytical approximations reported here and 



the numerical results reported in Capitan et al. (2009 ). 

For every R we can always find a community which is un- 
invadable at all its levels t < L. If 7? is such that ( (37] i is not 




Figure 5: Number of communities in the recurrent sets obtained with the analyt- 
ical approximations {N, showed with crosses) and with a numerical integration 
of the population dynamics (Wo, circles). The inset contains the absolute dif- 
ference \N - N(,\. The global picture is the same as that found in , Capitan et al.| 
j2009j, although differences of a few tens arise in some cases. 



verified, then the invader at level L + 1 initially decreases and 
goes extinct. This explains the intervals of R where only one ab- 
sorbent state is found. However, if ( (3T] i holds (with our choice 
of parameters this happens when si > 4), there is an initial time 
interval where the population of the invader is above the thresh- 
old. This can cause the extinction of lower level species, and 
generate recurrent sets with more than one community. 

Our analytical approximations thus provide results very close 
to those obtained numerically. Besides its being more efficient 
(the whole assembly needs not be generated), this method also 
allows to predict what would happen for values of R larger than 
1700, which are computationally prohibitive for the numerical 
method. With our bounds ( |24] l and ( (37] i we can estimate the 
next interval of R where more than one community in the end 
state will appear, namely R e [3844, 5 114]. That is out of reach 
of the numerical method, because the number of com munities 
in the whole assembly graph grows as fast asN ~ e""^ { Capitan 
[etaLilSoTOI . 



Two observations are on purpose. First, there are small inter- 
vals of R where the graph constructed starting from the empty 
community has L levels but there are viable communities with 
L + I levels which cannot be assembled starting from [this 
phenomenon is analogous to the existence of unreachable per- 
sistent communities showed in [Warren et al.| ( |2003] l]. We ob- 
serve this for R = 460, 465, 1615, 1620 and 1625 (see Table[TJ. 
We have checked that even in these cases the recurrent state is 
exactly recovered using the analytical approximations. 

Secondly, we can observe from Figures[5]and[6]that there are 
small regions where recurrent sets with more than one commu- 
nity are found out of the intervals predicted in Table [T] (around 
R 200 for L = 3 and R ^ 620 for L = 4). For those values, a 
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Figure 6: Number N of communities in tlie recurrent sets obtained with the 
analytical approximations (crosses in Figure [5}. The global t rend is the same as 
found in our previous work [see Figure 3. in |Capitan et al.| |2()09 1]. The inset 
shows the relative difference in the prediction of the number of communities in 
the end states. Note that the discrepancies occur in a region where this number 
is small. This explains the relatively large eiTor found in some cases. 



single absorbent community should be found. However, condi- 
tion ( [34] l for an invader at level L to initially grow and become 
extinct at equilibrium renders si < 2 for our choice of p. We 
have checked that this condition is satisfied by all these small 
recurrent sets, thus explaining their appearance. 

We have to assess the accuracy of the transitions predicted 
in the graph of our recurrent sets. Note that a slight difference 
in the ordering of extinctions can change the final community 
after the invasion and this may change the observed graph and 
therefore the asymptotic probability distribution of the associ- 
ated Markov chain. In order to check the transition matrices we 
obtain, we have calculated two averages. In Figure |7] we show 
the variation of the average number of species in the recurrent 
sets as a function ofR. The behavior is almost indistinguishable 
from that found in Capitan et al. ( 2009[ l (the inset of Figure |7] 
shows that the relative error is small). 

We have also checked that the number of extinctions pre- 
dicted with our approximations follows the same distribution 
than the one calculated numerically. To this purpose we de- 
fine the magnitude of an avalanche of extinctions as the relative 
variation m - AS /S of the total number of species in a com- 
munity after an invasion. In Figure [8] we show the cumulative 
histogram for the distribution of these magnitudes. We can see 
that the deviations between both distributions are small. Further 
statistical results will be discussed in the second paper of this 
suite. 



7. Conclusions 

In this paper we have presented a general model of trophic- 
level structured food-web, where interactions between species 




Figure 7: Average number of species S^y in the end states calculated analyti- 
cally vs. R. In the inset we show the relative error between S and its coiTe- 
sponding average So for each graph calculated numerically. 



are either feeding or competing. For the sake of simplicity, 
feeding only takes place between contiguous levels. The pop- 
ulation dynamics is modeled through Lotka-Volterra equations, 
and a proof is given that a wide class of these models has a glob- 
ally stable interior equilibrium. We have introduced this model 
as an appropriate general framework to study the process of 
successional invasions. In the invasion process, we consider a 
mean-field version, in which species in the same level are troph- 
ically equivalent and only intra- and interspecific competition is 
distinguished. This species symmetry assumption has allowed 
us to obtain analytical results, some of them exact and some 
other approximate. Among them we have provided estimations 
for the maximum number of species allowed per level, the max- 
imum number of levels for a given value of the resource satu- 
ration, and certain analytical approximations of the dependence 
of the equilibrium abundances on the occupancies of each level. 
We have combined these results with some criteria for the ac- 
ceptance of an invader in our model communities, and with the 
help of some global approximations of the invaded dynamics 
we have been able to obtain, with high accuracy, the sequence 
of extinctions occurring after an invasion. With this procedure 
we have reproduced the same results that we found in a pre- 
vious work ( Capitan et al.| [2009 ), this time without resorting 
to an integration of the Lotka-Volterra equations and without 
constructing the whole assembly graph. Among other things 
this brings the opportunity of exploring the model for resources 
which would otherwise be computationally prohibitive to ob- 
tain. 



Although the main results of this model are discussed at 
length in the second paper of this suite (Capitan et al. 2010| l, 
we have provided here a few of them which illustrate the global 
assembly process and some of its main features. For instance, 
we had reported already in Capitan et al. ( 2009| l that, upon in- 
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Figure 8: Cumulative probability function Il(m) for the distribution of the mag- 
nitude m of avalanches of extinctions. The distributions follow an exponential 
behavior. Crosses represent the results for our approximated transition ma- 
trix. The number of recurrent states coincide for the analytical and numerical 
method. The inset shows a case where the number of communities is underesti- 
mated. This explains the absence of several points in the distribution estimated 
analytically. The agreement is rather good even in this case. 



creasing the resource saturation R, the number of levels, L, that 
the system is able to sustain increases discontinuously. We pro- 
vide here an estimate of the values of R at which this occurs, 
and show that this values grow essentially as ~ ij+ly )^- Un- 
der the assumption that populations are close to the extinction 
level, we have shown that equilibrium communities are pyra- 
mids — again in agreement with the results obtained in Capitan 



et al. ( 2009 1. Close to the onset of appearance of a new level, 
the number of communities in the end state increases. We have 
identified that the requirement for this to happen is that the pop- 
ulation of a top predator invading the community initially grows 
only to go eventually extinct. From this knowledge we can es- 
timate the value of R at which the end state starts to have more 
than just one community. 

We have tested the approximations we have made by calcu- 
lating some observables. Among them we report on the average 
species richness as a function of R, as well as the distribution of 
the avalanche of extinctions produced by an invasion. In both 
cases the agreement is very good. In the latter case, it is worth 
mentioning that this distribution of avalanches decays exponen- 
tially with the avalanche size, meaning that there is a character- 
istic size of the avalanches. This size roughly grows with the 
species richness of the community, as one could expect. In any 
case, avalanches never get even close to destroy the community. 

We also propose in this paper an analytical approximation to 
the dynamics of a community invaded by a top predator. This 
approximation has been built matching the initial behavior of 
the solution (derived from the initial condition) and the asymp- 
totic decay expected close to the equilibrium. We have found a a 
rather good agreement with the solutions obtained by a numeri- 



cal integration of the Lotka-Volterra equations, and has allowed 
us to correctly predict (in most of the cases) the order of extinc- 
tions eventually caused by the invasion of a top predator. These 
approximations have been applied to reproduce the assembly 
graphs for the recurrent sets, showing small discrepancies only 
for certain values of R. This provides an alternative method to 
analyze the system for other sets of parameter values, with a 
negligible computational cost compared to the construction of 
the whole assembly graph. 

Our assembly model is based on several assumptions regard- 
ing the invasion process. Two of the most important ones are 
that newcomers invade at low population and the average time 
between invasions is large compared to the time for the com- 
munities to reach the equilibrium. If the invasion rate is too 
high ( Fukami 2004 Bastolla et al.' '2005 b| or if th e invasion 
is not produced by rare species ( Hewitt and Huxelj |2002[ ), the 
assembly process — and hence the resulting end states — can be 
drastically altered. The reason is that communities that are not 
accessible from the equilibrium state may be so from a tran- 
sient or if there is a massive invasion. This changes the assem- 
bly graph in ways that we can neither predict nor even check, 
because these processes are out of reach of our model. For in- 
stance, considering invading transients, one of the strong sim- 
plifications we make use of is that of starting always from a 
well-defined initial condition, namely the equilibrium state. If 
the system can be invaded at any moment during a transient 
there are infinitely many initial conditions to start off from, 
something we cannot implement. So what happens if any of 
those two hypotheses is violated remains an open question. 
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Appendix A. Derivation of the reduced dynamical system 

We will show in this appendix that our dynamical system 
«^ - q^n^^, with the linear response function ([12]), can be re- 
duced to the form ( fTSj l when all the initial species abundances 
at a certain level are equal. The crucial point for this to be true 
is the relation ([14]). 

This result can be formulated in a simple way. Consider the 
two-dimensional autonomous system 



(A. I) 



with the initial condition x{Q) - y{0) and which satisfies 
f{x,y) - g{y,x). We are going to show that the Taylor ex- 
pansions centered at f = of x and y are identical. In principle, 
both expansions will have certain radii of convergence. Let t be 
lower than the minimum of these radii. Then we just need to 
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show that all the derivatives at f = coincide. But this follows 
by induction. 

The first derivatives are shown to be equal easily. Let us 
assume that jcW(O) = yW(0) for all ;t = Then the 

(n + l)-th derivative is 



d"f 



x^j\Q)y^"-j\0). (A.2) 



dxjdy"-J 

But, since f{x, y) = g(y, x), this is equivalent to write 



^ W / dyJdx"-J 



and, relabelling the sum index. 



(A.3) 



t=0 



dxidy"-i 



x^j\<d)y'-"-J\Q), (A.4) 



f=0 



which is equal to /"+'\0). 

Therefore we have shown that the Taylor expansions of x(f) 
and y{t) coincide. This means that x{t) - y{t) within the radius 
of convergence of the series. For larger times, we can apply 
the same argument by analytic continuation (we choose some 
fo in the interval or convergence as the centering point for a 
new Taylor expansion, and repeat the argument). Hence we 
conclude that x{t) - y{t) for all t. 

Note that the same considerations apply to our system ( [T2] l, 
so we can reduce considerably the complexity of the system and 
solve ([TSj instead. 

Appendix B. Analytical approximation to the equilibrium 
densities 

This appendix is devoted to solve the linear system for the 
equilibrium densities ( |25] l. The solution of this system can be 
obtained through Cramer's rule as 



(B.l) 



for certain determinants H/./ and A/.. Our approximation is 
based in some recurrent equations that can be obtained for these 
determinants. 

Let us start with the (L+ l)x(L+ 1) determinant 



(B.2) 



where d[ = p + Hence the densities depend on {s(}^^^ 
only through the inverse of all the possible products , ■ ■ ■ , 
for some combination (11,12, ■■■ Jk) of k elements of the set 
{1,2, ... ,L}. In the recurrent sets we get the highest occupancy 
of species in each level allowed by the resource according to 
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((2T|-((23ll, so we expect that a rather good approximation for 
the equilibrium densities amounts to neglecting orders higher 
than 1/i. Hence 



(B.3) 



where 
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(B.4) 
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has dimension (L+l)x(L+l) and Bl( is the determinant 
obtained by substituting the ^-th column of Di by the column 
vector ug whose components are , = 5f , (for / = 0, 1, . . . L). 
The determinant De satisfies the recursion 



(B.5) 



where £ - 1,2, . . .L, Do — -\ and Di - p + y+j-. This relation 
can be easily solved using a generating function. On the other 
hand, it is easy to see that Bi c - Df_i£'£_f_i, with E{ the {{ + 
1) X (^ + 1) determinant 



-p 


-y- 


■ 


■ 


y+ 


-p 


-y- ■ 


■ 





y+ 


-p ■ 


■ 








■ 


■ -p 



(B.6) 



which also satisfies recursion ( |B.5| l with Eq - -p and Ey 

p^ + y+y-- 

The generating function that results from (|B.5|l is 



(=0 



A) + (£>i +pDo)z 

y+y-z^ - pz- I 



(B.7) 



(B.8) 



with ae given by ^TEf . Then the following compact expressions 
result 



and after the series expansion we get 

D( = i-y-Y-' [(Di +pDo)ae-i - y-Eoae] , 



Of = (-irV-[flf + r+flf-i], 



(B.9) 
(B.IO) 



The explicit expression for H/,/ is obtained from A^ sub- 
stituting its f-th column by the (L + 1) x 1 column vector 
(-R, a,.. ., a)^. We can expand it up to leading order in powers 
of 1 /i to get 

^ / 1 \ 

EL,e^n,f-(l-p)Y,^+Oi-^\. (B.ll) 
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where 
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and 2^ ^ is the determinant that results when we substitute the 
j-th column of Tl ( by uj {j + t). 

Expanding T/,/ along its first row we get 

Tut = -aA^,i + aj^j-AL-u-i + {-\i^^Rj\E^-t-x, (B.13) 
where we define the new / x ; determinants A,j as 
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that satisfy the recurrence equation 

A„j = -pA„-ij + y+y-A„-2j + y"S'Ej-2, 



for j - 1,2, ... ,n - I (with the boundary conditions Ajj+i - 
and Ajfi = 0), and 



An,n - -y+A„-l^„-\ + En-2- 



(B.16) 



These relations can be explicitly solved. On the one hand, by 
definition Al l = 1, which amounts to choosing £-1 = 1 for this 
to be compatible with ( |B.16| l. Moreover, making use again of a 
generating function, the solution of ( |B.16[ ) is 



r- - r+ + p 



p ^ r+ -f 

— a^-i + 

y- y- 



p ^ r+ / 74 

— flj_2 + — - — 

y- \y- 



(B.17) 



for i >2. On the other hand, the explicit solution of ( |B.15| l is 

Aj+k.j = (-l)*r-'"'at+iAjj 

+ \{-\f^\y-a, - r+fli-i) + r- , 

y- - 7+ + p L J 

for k>\. Therefore equations ( |B.17| i and ( |B.18| ), together with 
( |B.13[ ), provide an explicit solution for the determinants T^ f. 

Fortunately, f can be written in terms of the previous de- 
terminants since Qi i is a block-diagonal determinant with two 
blocks that satisfies 

e{, = D,_,A^,,,_,., for k<j, (B.19) 

Q[i^EL-j-iTj^U^ for k>j. (B.20) 

This completes the analytical approximation of the equilibrium 
densities of our dynamical model. We have derived explicit ex- 
pressions for all the terms involved in ( |B.l[ l, ( |B.3[ ) and ( |B.l 1[ ) up 



to leading order in 1 1 s. Moreover, note that the same technique 
applied to find this approximation can be extended to obtain the 
exact dependence on {^fl^^j of the abundances. Higher-order 
terms in powers of 1 /i introduce in the corresponding deter- 
minants several column vectors of the type of ui making each 
determinant to be block-diagonal involving Di, E(, Ai j or T^^, 
so that the general solution contains in each term a product of 
a certain combination of these determinants. This explicit ex- 
pression can in fact be written, but it is too cumbersome. The 
approximations here obtained are both sufficiently simple and 
accurate enough to capture the behavior of population densities 
in the communities of the recurrent sets. 

Appendix C. Technical details of the global approxima- 
tions to the dynamics 

In this appendix we will describe the calculation of the unde- 
termined parameters of our ansatz (|35|)-((36l) for the dynamics 
of system invaded by a top predator We impose that the initial 
condition and the first k derivatives at f = match the exact 
values, which can be readily calculated. Indeed, our system has 
the form i,- - -o-x,- -i- xifj{x), where fiix) - 2; bijXj is a linear 
function. Therefore we can recursively calculate the s + \ initial 
derivative as 



(B.15) x<"'>(0) = -a4"(0) + 



(0)/i(xW(0)). (C.l) 



For a real eigenvalue {cj - 0), we choose C(f) [see ([36|] to 
be a polynomial of degree k-2, and for a complex one (oj + 0) 
we choose degree A; - 3, in order to compensate for the extra 
undetermined coefficient in the long-term behavior in this case. 
Equating the approximate solution to the initial condition and 
the first k-\ derivatives of our ansatz to the exact values leads 
to a linear system for the undetermined coefficients. The equa- 
tion for the ^-th derivative yields a polynomial equation for ^, 
namely 



when w = 0, where 



|](^-2)///->-2 = (^Hc.V^', (c. 



2) 



Hj = + o?-)n^'\Q) + 2/in*-'+"(0) + n*-'+2>(0) (C.3) 

and n^'^ stands for the j-th derivative of n^^', which can be cal- 
culated exactly using ( |C.1| ). For u) - ^ Eq. ( |C.2| ) gets replaced 
by 

Y\^~-^\ + n^^^'m] t^-' - ^P^"'- (C.4) 

Afterwards, we just need to calculate the coefficients cj and do 
(and t/i , if dLt 0) by solving the linear system that they satisfy. 

Once we have the approximate time behavior for n^^' we cal- 
culate analytically the remaining populations by direct sub- 
stitution into the system ^T5\ , taking advantage of the recursive 
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form of these equations, once n^^' is known. Notice that, since 
we have to calculate successive derivatives in order to get any 
lower population, the accuracy of «^^' at short times degrades 
as we calculate lower level populations. Fortunately the model 
produces communities with a small number of trophic levels 



(Capitan et al. 2009 1. The choice ^ = 5 seems to be enough 
to account for the dynamics of any community of up to L = 4 
levels invaded by a top predator (see Figures |2] and [3]). For the 
description of the dynamics of communities with a higher num- 
ber of levels we would need to choose polynomials of higher 
degree in our ansatz. 

A final caveat needs to be made with respect to the calcu- 
lation of ^. We need it to be positive, otherwise ([36| would 
be meaningless. Among all the roots of ( |C.2| l we choose the 
largest, positive, real solution, so that any possible initial oscil- 
lation of the polynomial C(f) is damped by the exponential. In 
the majority of the dynamics that we have approximated (see 
Section |6]l, we are able to find a positive solution for ^. How- 
ever, in some cases there is no positive solution. In those cases 
we just minimize the difference between the exact k-th deriva- 
tive and the approximate one at f = 0. This also produces an ac- 
ceptable solution. In all minimization procedures that we have 
run, a positive exponent ^ is always found. 
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